#include <stdint.h>
#include <math.h>
#include <assert.h>
#include <limits.h>

#include "myexp.h"
#include "myexp_impl.h"
#include "myexp.inl"

static double specialcase(double x)
{
    uint32_t abstop = top12(x) & 0x7ff;
    if (unlikely(abstop >= 10 + DBL_EXP_BASE)) // x <=-1024(2^10) or x >= 1024.0
    {
        if (isnan(x))
        {
            return NAN; // return nan
        }
        if (x > 0)
            return __builtin_inff();
        else
            return 0.0;
    }

    double_t z, z0, r, r2, r4, body, tail;
    double scale = 1.0;
    z = x * INV_LN2; // z = x / ln(2)
    union
    {
        double f;
        int32_t i[2];
    } zi = {z + 0x1.8p45}; // zi.i[0] = (round)(z *128)

    int k = (zi.i[0] >> 7); // here, -1477 <= k <= 1477
    int i = zi.i[0] - k * 128;

     // here, -1477 <= k <= 1477,0 <= idx < 128
    if (k > 1023)
    {
        scale = 0x1p1023;
        k -= 1023;
    }
    //The following code require k-61 >= -1022, so k >= -961
    else if (k < -961)
    {
        scale = 0x1p-961;
        k -= (-961);
    }

    //-1477-(-961)=-516,1477-(1023) =454
    // so -961 <= k<= 1023
    // z0是z的近似值，z0具有这样的形式，z0 = k + i/128.0, -961<k<1024, 0 <= i < 128
    z0 = zi.f - 0x1.8p45;
    r = x - (z0 * Ln2hiN) - (z0 * Ln2loN); // 将ln(2)拆分为２个数，减少误差, r= x - (i+k/128)*ln(2)
    r2 = r * r;
    r4 = r2 * r2;

    int64_t i_body = exp_data_61[i];                    // i_body =power(2.0,i/128) * 2^61
    tail = r2 * (C2 + r * C3) + r4 * (C4 + r * C5) + r; // tail= exp(r)-1.0

    // 转化浮点数为定点数格式，i_tail包括61位小数
    int64_t i_tail = double_2_fixpoint(tail);

    // 在计算res时，使用定点数方法来计算，使得中间计算过程的精度超过double浮点数的精度
    // res= body*exp(r) = body*(1+tail)= body+body*tail
    int64_t res = i_body + fixpoint_mul(i_body, i_tail);

    // 当 k 满足条件 -1022 <= k-61 < 1024 时, pow2(k-61)是规格化的双精度浮点数
    // -1022+61 = -961 <= k < 1024+61
    return (double)res * pow2(k - 61) * scale;
}
/* 求exp(x)的核心算法
 基本原理，
  令 e^x　=2^z            (1)
  两边取对数得 x = z*ln  　 (2）
  =>　       z = x/ln(2)
   则　e^x = 2^z= 2^(x/ln(2))
   我们取z的近似值为z0, z0具有这样的格式, z0=k + i/128, k是整数,i是非负整数，0<=i<128．
   z0满足　z-1/256 <= z0 <= z+1/256, 即abs(z-z0)<＝1/256
   z0的计算过程为：
   　1. t = round(z*128.00), 将z放大128倍，并四舍五入取整，这里，round(z*128.00)是利用浮点数的特点实现的，稍后解释
   　2. k = (t>>7), t>>7等价于t除以128向负无穷取整，
      当t是一个非负整数时，t除以128的商(记作t/128)等于t>>7,
      当t就一个负数时，t/128不等于(t>>7), 如-129/128=-1,而 -129>>7=-2，使用移位运算，保证了以下的i总是非负整数
    3. i  = t-k*128,
    4. z0 = k + i/128
   下面的代码，对z0的计算是利用浮点数的特点实现的，稍后解释
　这个版本中，为了保证精度，部分代码是采用定点数来计算的．我们采用64位带符号整数(补码格式）表示一个定点数，
　包括1位符号位，2位整数位，61位小数位构成．当最高比特位为0时，这个定点数是一个正数或0，当最高比特位为１，这个定点数是一个负数．
  次高位的权重为2，最低位的权重是2^-61,这种表示方式，可表示区间的数的范围是开区间(-4,4)
*/

double exp_v3(double x)
{
    uint64_t abs_u64 = as_uint64(x) & LLONG_MAX;
    if (unlikely(abs_u64 > LN2_X_961_AS_UINT64)) // x > 961*ln(2) or x < -961*ln(2)
        return specialcase(x);

    double z = x * INV_LN2; // z = x / ln(2)
    union
    {
        double f;
        int32_t i[2];
    } zi = {z + 0x1.8p45}; // zi.i[0] = (round)(z *128)

    // here -961 <= k <= 961, 0 <= i < 128
    int k = (zi.i[0] >> 7);     
    int i = zi.i[0] - k * 128;

    // z0是z的近似值，z0具有这样的形式，z0 = k + i/128.0, 
    double z0 = zi.f - 0x1.8p45;
    double r = x - (z0 * Ln2hiN) - (z0 * Ln2loN); // 将ln(2)拆分为２个数，减少误差, r= x - (i+k/128)*ln(2)
    double r2 = r * r;
    double r4 = r2 * r2;

    int64_t i_body = exp_data_61[i];                    // i_body =power(2.0,i/128) * 2^61
    double tail = r4 * (C4 + r * C5) + r2 * (C2 + r * C3) +  r; // tail= exp(r)-1.0

    // 转化浮点数为定点数格式，i_tail包括61位小数
    int64_t i_tail = double_2_fixpoint(tail);

    // 在计算res时，使用定点数方法来计算，使得中间计算过程的精度超过double浮点数的精度
    // res= body*exp(r) = body*(1+tail)= body+body*tail
    int64_t res = i_body + fixpoint_mul(i_body, i_tail);

    // 当 k 满足条件 -1022 <= k-61 < 1024 时, pow2(k-61)是规格化的双精度浮点数
    // -1022+61 = -961 <= k < 1024+61
    return (double)res * pow2(k - 61);
}
